R Notebook

library(tidyverse)

Loading tidyverse: ggplot2
Loading tidyverse: tibble
Loading tidyverse: tidyr
Loading tidyverse: readr
Loading tidyverse: purrr
Loading tidyverse: dplyr
Conflicts with tidy packages -----------------------------------------------------------------
filter(): dplyr, stats
lag():    dplyr, stats

library(vcfR)

   *****       ***   vcfR   ***       *****
   This is vcfR 1.5.0 
     browseVignettes('vcfR') # Documentation
     citation('vcfR') # Citation
   *****       *****      *****       *****

format.names <- str_split(data$FORMAT[[1]],pattern=":")[[1]]
data2 <- data %>%
  select(-INFO,-ID,-FILTER, -FORMAT) %>%
  mutate(S95=str_replace_all(S95,"([^0-9]|$)\\.","\\1NA")) %>% # convert missing data "." into NA
  mutate(S96=str_replace_all(S96,"([^0-9]|$)\\.","\\1NA")) %>%
  separate(S95,sep=":",into = str_c("S95_",format.names), convert = TRUE) %>%
  separate(S96,sep=":",into=str_c("S96_",format.names), convert = TRUE) %>%
  filter(!(S95_GT=="1/1" & S96_GT=="1/1")) %>%
  filter(!(S95_GT=="0/0" & S96_GT=="0/0")) %>%
  select(-ends_with("DPR"))
data2

calculate allele frequencies and polarize based on S95 for each SNP

data2 <- data2 %>%
  mutate(S95_alt_AF=S95_AO/S95_DP,
         S96_alt_AF=S96_AO/S96_DP,
         S95_alt_major=S95_alt_AF > 0.5)
data2 <- data2 %>%
  mutate(S95_AF_polar=ifelse(S95_alt_major,S95_alt_AF,1-S95_alt_AF),
         S96_AF_polar=ifelse(S95_alt_major,S96_alt_AF,1-S96_alt_AF))
data2         

rolling mean

plot!

plA.data <- pl.data2 %>% filter(str_detect(CHROM,"chrA"),!str_detect(CHROM,"random")) 

plA <- ggplot(plA.data,aes(x=POS,y=AF_100,color=population,shape=population)) +
  facet_wrap(~CHROM, scales = "free_x",ncol=2) +
  #geom_point(alpha=.1) +
  geom_line()
plA

plot!

plC.data <- pl.data2 %>% filter(str_detect(CHROM,"chrC"),!str_detect(CHROM,"random")) 

plC <- ggplot(plC.data,aes(x=POS,y=AF_100,color=population,shape=population)) +
  facet_wrap(~CHROM, scales = "free_x",ncol=2) +
  #geom_point(alpha=.1) +
  geom_line()
plC

plot one per chromosome

chrom.names <- unique(pl.data2$CHROM) %>% str_subset("[0-9]$")
chrom.names

 [1] "chrA01" "chrA02" "chrA03" "chrA04" "chrA05" "chrA06" "chrA07" "chrA08" "chrA09" "chrA10" "chrC01" "chrC02" "chrC03" "chrC04"
[15] "chrC05" "chrC06" "chrC07" "chrC08" "chrC09"

for (chr in chrom.names) {
  pl.data2 %>% filter(str_detect(CHROM,chr)) %>%
    ggplot(aes(x=POS,y=AF_100,color=population)) +
    geom_line() +
    ggtitle(chr) +
    ylim(0,1)
  ggsave(str_c(chr,".png"),height = 8, width = 12)
}

absolute difference in allele frequency

data3 <- data2 %>%
  group_by(CHROM) %>%
  mutate(abs_delta_af = abs(S95_AO/S95_DP - S96_AO/S96_DP),
         abs_delta_af_100=roll_mean(abs_delta_af,n=100,fill=NaN),
         abs_delta_af_500=roll_mean(abs_delta_af,n=500,fill=NaN)
         )
data3

plot one per chromosome

for (chr in chrom.names) {
  data3 %>% filter(str_detect(CHROM,chr)) %>%
    ggplot(aes(x=POS,y=abs_delta_af_100)) +
    geom_line() +
    ggtitle(chr) +
    ylim(0,1)
  ggsave(str_c("Delta_af",chr,".png"),height = 8, width = 12)
}

Try filtering for F2 good snps

F2snps <- read_tsv("Final_F2_SNP_Sites.tab")

Parsed with column specification:
cols(
  CHROM = col_character(),
  POS = col_integer()
)

data2.F2.filtered <- semi_join(data2,F2snps)

Joining, by = c("CHROM", "POS")

dim(data2)

[1] 1926088      26

dim(F2snps)

[1] 18226     2

dim(data2.F2.filtered)

[1] 1997   26

data3.F2 <- data2.F2.filtered %>%
  group_by(CHROM) %>%
  mutate(abs_delta_af = abs(S95_AO/S95_DP - S96_AO/S96_DP),
         abs_delta_af_10=roll_mean(abs_delta_af,n=10,fill=NaN)
         )
data3.F2

for (chr in chrom.names) {
  data3.F2 %>% filter(str_detect(CHROM,chr)) %>%
    ggplot(aes(x=POS,y=abs_delta_af_10)) +
    geom_line() +
    ggtitle(chr) +
    ylim(0,1)
  ggsave(str_c("Delta_af_F2_filter",chr,".png"),height = 8, width = 12)
}

filter on 505

snps.505 <- read_csv("../kiat/Breeding_Prediction_App/Data/505/505_Genotype_170K_SNPs.csv.gz") %>%
  select(CHROM,POS)

Duplicated column names deduplicated: 'ID_505_K200' => 'ID_505_K200_1' [21]Parsed with column specification:
cols(
  .default = col_integer(),
  MARKER = col_character(),
  CHROM = col_character(),
  REF = col_character(),
  ALT = col_character()
)
See spec(...) for full column specifications.

head(snps.505)

data2.505.filtered <- semi_join(data2,snps.505)

Joining, by = c("CHROM", "POS")

dim(data2)

[1] 1926088      26

dim(snps.505)

[1] 174397      2

dim(data2.505.filtered)

[1] 14116    26

data3.505 <- data2.505.filtered %>%
  group_by(CHROM) %>%
  mutate(abs_delta_af = abs(S95_AO/S95_DP - S96_AO/S96_DP),
         abs_delta_af_10=roll_mean(abs_delta_af,n=10,fill=NaN)
         )
data3.505

for (chr in chrom.names) {
  data3.505 %>% filter(str_detect(CHROM,chr)) %>%
    ggplot(aes(x=POS,y=abs_delta_af_10)) +
    geom_point(aes(y=abs_delta_af),shape=1) +
    geom_line(color="blue") +
    ggtitle(chr) +
    ylim(0,1)
  ggsave(str_c("Delta_af_505_filter",chr,".png"),height = 8, width = 12)
}